4DN DamID - Ki67 project - Actinomycin D analysis
Tom van Schaik
2022-05-06
Project
Mapping of Ki67 interactions with the genome and comparison with lamina interactions.
Introduction
Various analyses of cells treated with low doses of actinomycin D or DMSO to perturb nucleolar targeting of Ki67.
Method
NA
Set-up
Set the parameters and list the data.
# Load dependencies - without warnings / messages
library(tidyverse)
library(GenomicRanges)
library(rtracklayer)
library(ggplot2)
library(RColorBrewer)
library(GGally)
library(corrr)
library(caTools)
library(colorspace)
library(ggrastr)
# Prepare output
output_dir <- "ts220503_6_ki67_actinomycin"
dir.create(output_dir, showWarnings = FALSE)
# Load input
chromosomes <- c(paste0("chr", 1:22), "chrX")
input_dir <- "ts220503_1_data_gathering"
bin_size <- readRDS(file.path(input_dir, "bin_size.rds"))
centromeres <- readRDS(file.path(input_dir, "centromeres.rds"))
colors_set1 <- readRDS(file.path(input_dir, "colors_set1.rds"))
colors_set2 <- readRDS(file.path(input_dir, "colors_set2.rds"))
colors_set3 <- readRDS(file.path(input_dir, "colors_set3.rds"))
tib_padamid_replicates <- readRDS(file.path(input_dir,
"tib_padamid_replicates.rds"))
tib_padamid_combined <- readRDS(file.path(input_dir,
"tib_padamid_combined.rds"))
tib_padamid_combined_scaled <- readRDS(file.path(input_dir,
"tib_padamid_combined_scaled.rds"))
gr_padamid_replicates <- readRDS(file.path(input_dir,
"gr_padamid_replicates.rds"))
gr_padamid_combined <- readRDS(file.path(input_dir,
"gr_padamid_combined.rds"))
tib_hmm_replicates <- readRDS(file.path(input_dir, "tib_hmm_replicates.rds"))
tib_hmm_combined <- readRDS(file.path(input_dir, "tib_hmm_combined.rds"))
padamid_metadata_replicates <- readRDS(file.path(input_dir,
"padamid_metadata_replicates.rds"))
padamid_metadata_combined <- readRDS(file.path(input_dir,
"padamid_metadata_combined.rds"))
input_dir <- "ts220503_5_osmotic_shock"
tib_damid_combined <- readRDS(file.path(input_dir, "tib_damid_combined.rds"))
# Scale pA-DamID scores?
tib_padamid_combined <- tib_padamid_combined %>%
mutate_at(4:ncol(.), function(x) scale(x)[, 1]) %>%
filter(seqnames != "chrY")library(knitr)
opts_chunk$set(fig.width = 10, fig.height = 4, cache = T,
dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/"))
pdf.options(useDingbats = FALSE)# From Fede:
# ggpairs custom functions
corColor <- function(data, mapping, color = I("black"), sizeRange = c(1, 3), ...) {
x <- eval_data_col(data, mapping$x)
y <- eval_data_col(data, mapping$y)
r <- cor(x, y)
rt <- format(r, digits = 3)
tt <- as.character(rt)
cex <- max(sizeRange)
# helper function to calculate a useable size
percent_of_range <- function(percent, range) {
percent * diff(range) + min(range, na.rm = TRUE)
}
# plot correlation coefficient
p <- ggally_text(label = tt, mapping = aes(), xP = 0.5, yP = 0.5,
# size = I(percent_of_range(cex * abs(r), sizeRange)),
size = 6,
color = color, ...) +
theme(panel.grid.minor=element_blank(),
panel.grid.major=element_blank())
corColors <- RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")[2:6]
if (r <= boundaries[1]) {
corCol <- corColors[1]
} else if (r <= boundaries[2]) {
corCol <- corColors[2]
} else if (r < boundaries[3]) {
corCol <- corColors[3]
} else if (r < boundaries[4]) {
corCol <- corColors[4]
} else {
corCol <- corColors[5]
}
p <- p +
theme(panel.background = element_rect(fill = corCol))
return(p)
}
customScatter <- function (data, mapping)
{
p <- ggplot(data = data, mapping) +
geom_bin2d(bins = 100) +
geom_smooth(method = "lm", se = T, col = "red") +
scale_fill_gradient(low = "lightgrey", high = "black", name = "Count") +
theme_bw()
p
}
# Plots
PlotDataTracksLight <- function(input_tib, exp_names, centromeres,
color_groups, plot_chr = "chr1",
plot_start = 1, plot_end = 500e6,
smooth = 1, color_list = NULL,
fix_yaxis = F, aspect_ratio = NULL,
lighten_negative = T, raster = T,
ylimits = NULL) {
# Get the scores for the samples
tib_plot <- input_tib %>%
dplyr::select(seqnames, start, end,
all_of(exp_names))
if (smooth > 1) {
tib_plot <- tib_plot %>%
mutate_at(vars(contains("_")),
runmean, k = smooth)
}
# Filter for plotting window
tib_plot <- tib_plot %>%
filter(seqnames == plot_chr,
start >= plot_start,
end <= plot_end)
# Gather
tib_gather <- tib_plot %>%
gather(key, value, -seqnames, -start, -end) %>%
mutate(key = factor(key, levels = exp_names),
fill_column = color_groups[match(key, names(color_groups))],
fill_column = factor(fill_column, levels = unique(color_groups)))
# If desired, make negative values a lighter shade to improve visualization
if (lighten_negative) {
tib_gather <- tib_gather %>%
mutate(fill_column = interaction(fill_column,
value < 0))
}
# Centromeres
centromeres.plt <- as_tibble(centromeres) %>%
filter(seqnames == plot_chr) %>%
gather(key_centromeres, value, start, end)
# Plot
if (is.null(ylimits)) {
ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
}
fill_column <- NULL
plt <- tib_gather %>%
ggplot(aes(fill = fill_column))
# Set all counts tracks to the same limits
if (raster) {
plt <- plt +
rasterize(geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6,
ymin = 0, ymax = value)),
dpi = 300)
} else {
plt <- plt +
geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6,
ymin = 0, ymax = value))
}
plt <- plt +
geom_hline(yintercept = 0, size = 0.5) +
geom_line(data = centromeres.plt,
aes(x = value / 1e6, y = 0),
color = "black", size = 2) +
facet_grid(key ~ ., scales = "free_y") +
xlab(paste0(plot_chr, " (Mb)")) +
ylab("Score") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0.025, 0.025)) +
theme_classic()
if (! is.null(color_list)) {
if (lighten_negative) {
color_list <- c(color_list,
lighten(color_list, amount = 0.5))
}
colors <- levels(tib_gather$fill_column)
color_list <- color_list[1:length(colors)]
names(color_list) <- colors
#colors <- recode(colors, !!!color_list)
plt <- plt +
scale_fill_manual(values = color_list, guide = "none")
} else {
plt <- plt +
scale_fill_brewer(palette = "Set1", guide = "none")
}
if (fix_yaxis) {
plt <- plt +
coord_cartesian(xlim = c(max(c(plot_start,
min(tib_gather$start))) / 1e6,
min(c(plot_end,
max(tib_gather$end))) / 1e6),
ylim = ylimits)
} else {
plt <- plt +
coord_cartesian(xlim = c(max(c(plot_start,
min(tib_gather$start))) / 1e6,
min(c(plot_end,
max(tib_gather$end))) / 1e6))
}
if (! is.null(aspect_ratio)) {
plt <- plt +
theme(aspect.ratio = aspect_ratio)
}
plot(plt)
}
PlotScatterBinned <- function(tib, n1, n2, color_by = NULL, identity = F,
n_min = 10, ylimits_col = c(-2.4, 2.4),
count_range = c(0, 400), smooth = 1,
remove_outliers = F) {
# Get tibble
tib_process <- tib %>%
dplyr::select(seqnames, all_of(n1), all_of(n2)) %>%
rename_all(~ c("seqnames", "n1", "n2"))
# Calculate pearson correlation without any smoothing
tib_cor <- round(cor(tib_process$n1, tib_process$n2, use = "complete.obs",
method = "pearson"), 2)
if (smooth > 1) {
tib_process <- tib_process %>%
group_by(seqnames) %>%
mutate(n1 = runmean(n1, k = smooth),
n2 = runmean(n2, k = smooth))
}
if (! is.null(color_by)) {
tib_process <- tib_process %>%
add_column(color = color_by)
}
tib_process <- tib_process %>%
drop_na()
# Change color range
if (! is.null(color_by)) {
limits_color <- quantile(tib_process$color, c(0.001, 0.999), na.rm = T)
tib_process$color[tib_process$color < limits_color[1]] <- limits_color[1]
tib_process$color[tib_process$color > limits_color[2]] <- limits_color[2]
}
# Replace outliers
if (remove_outliers) {
tib_process <- tib_process %>%
mutate(n1 = case_when(n1 > quantile(n1, 0.999) ~ quantile(n1, 0.999),
n1 < quantile(n1, 0.001) ~ quantile(n1, 0.001),
T ~ n1),
n2 = case_when(n2 > quantile(n2, 0.999) ~ quantile(n2, 0.999),
n2 < quantile(n2, 0.001) ~ quantile(n2, 0.001),
T ~ n2))
}
# Metrics
n1_min = min(tib_process$n1) - 0.001
n1_max = max(tib_process$n1) + 0.001
n1_binsize <- (n1_max - n1_min) / 49
n2_min = min(tib_process$n2) - 0.001
n2_max = max(tib_process$n2) + 0.001
n2_binsize <- (n2_max - n2_min) / 49
tib_summary <- tib_process %>%
mutate(n1_cut = cut(n1,
seq(n1_min, n1_max, length.out = 50)),
n2_cut = cut(n2,
seq(n2_min, n2_max, length.out = 50))) %>%
mutate(n1_bin = as.numeric(as.factor(n1_cut)),
n2_bin = as.numeric(as.factor(n2_cut))) %>%
mutate(n1_bin = n1_min - n1_binsize/2 + n1_bin * n1_binsize,
n2_bin = n2_min - n2_binsize/2 + n2_bin * n2_binsize) %>%
group_by(n1_bin, n2_bin)
# Plot
if (! is.null(color_by)) {
tib_summary <- tib_summary %>%
dplyr::summarise(n = n(),
mark = mean(color)) %>%
ungroup() %>%
filter(n >= n_min)
plt <- tib_summary %>%
ggplot(aes(x = n1_bin, y = n2_bin)) +
geom_tile(aes(fill = mark)) +
scale_fill_gradient2(low = "blue", mid = "grey", high = "red",
midpoint = 0, limits = ylimits_col,
na.value = "green")
} else {
tib_summary <- tib_summary %>%
dplyr::summarise(n = n()) %>%
ungroup() %>%
filter(n >= n_min)
plt <- tib_summary %>%
ggplot(aes(x = n1_bin, y = n2_bin)) +
geom_tile(aes(fill = n)) +
scale_fill_gradient(low = "lightgrey", high = "black", name = "Count",
limits = count_range, na.value = "green")
}
plt <- plt +
geom_hline(yintercept = 0, linetype = "dashed", col = "black") +
geom_vline(xintercept = 0, linetype = "dashed", col = "black") +
xlab(n1) +
ylab(n2) +
ggtitle(paste0("Pearson: ", tib_cor)) +
theme_bw() +
theme(aspect.ratio = 1)
if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0,
col = "black", linetype = "dashed")
plot(plt)
}1. Correlations between time points
I will perform these analyses only on the combined data, as the data is only shallowly sequenced.
# Gather metadata
padamid_metadata_replicates_act <- padamid_metadata_replicates %>%
filter(experiment %in% c("actinomycin") &
target == "Ki67")
padamid_metadata_combined_act <- padamid_metadata_combined %>%
filter(experiment %in% c("actinomycin") &
target == "Ki67")
# Use GGally to make correlation plots
boundaries <- seq(from = 0.1, to = 0.7, length.out = 4)
# Replicates
tib <- tib_padamid_replicates %>%
drop_na() %>%
dplyr::select(padamid_metadata_replicates_act$sample)
plt <- ggpairs(tib,
upper = list(continuous = corColor),
lower = list(continuous = customScatter),
diag = list(continuous = function(data, mapping, ...) {
ggally_densityDiag(data = data, mapping = mapping, alpha = 0.3, fill = "red") +
theme_bw()})) +
ggtitle("Correlating all data") +
xlab("") +
ylab("")
print(plt)
# Combined
tib <- tib_padamid_combined %>%
drop_na() %>%
dplyr::select(padamid_metadata_combined_act$sample)
plt <- ggpairs(tib,
upper = list(continuous = corColor),
lower = list(continuous = customScatter),
diag = list(continuous = function(data, mapping, ...) {
ggally_densityDiag(data = data, mapping = mapping, alpha = 0.3, fill = "red") +
theme_bw()})) +
ggtitle("Correlating all data") +
xlab("") +
ylab("")
print(plt)
PlotDataTracksLightFaceted <- function(input_tib, exp_names, centromeres,
color_groups, facet_groups,
plot_chr = "chr1",
plot_start = 1, plot_end = 500e6,
smooth = 1, color_list = NULL,
fix_yaxis = F) {
# Exp names is a vector of sample names
exp_search <- paste(exp_names, collapse = "|")
# Get the scores for the samples
tib_plot <- input_tib %>%
dplyr::select(seqnames, start, end,
matches(exp_search))
if (smooth > 1) {
tib_plot <- tib_plot %>%
mutate_at(vars(contains("_")),
runmean, k = smooth)
}
# Filter for plotting window
tib_plot <- tib_plot %>%
filter(seqnames == plot_chr,
start >= plot_start,
end <= plot_end)
# Gather
tib_gather <- tib_plot %>%
gather(key, value, -seqnames, -start, -end) %>%
mutate(key = factor(key, levels = exp_names),
fill_column = color_groups[match(key, names(color_groups))],
fill_column = factor(fill_column, levels = unique(color_groups)),
facet_column = facet_groups[match(key, names(facet_groups))],
facet_column = factor(facet_column, levels = unique(facet_groups)))
# Centromeres
centromeres.plt <- as_tibble(centromeres) %>%
filter(seqnames == plot_chr) %>%
gather(key_centromeres, value, start, end)
# Plot
ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
fill_column <- NULL
plt <- tib_gather %>%
ggplot(aes(fill = fill_column))
# Set all counts tracks to the same limits
plt <- plt +
#geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6,
# ymin = 0, ymax = value)) +
geom_line(aes(x = (start+end)/2e6, y = value,
col = fill_column)) +
geom_line(data = centromeres.plt,
aes(x = value / 1e6, y = 0),
color = "black", size = 2) +
geom_hline(yintercept = 0, size = 0.5) +
facet_grid(facet_column ~ ., scales = "free_y") +
xlab(paste0(plot_chr, " (Mb)")) +
ylab("Score") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0.025, 0.025)) +
theme_classic()
if (! is.null(color_list)) {
colors <- levels(tib_gather$fill_column)
color_list <- color_list[1:length(colors)]
names(color_list) <- colors
#colors <- recode(colors, !!!color_list)
plt <- plt +
scale_fill_manual(values = color_list, guide = F) +
scale_color_manual(values = color_list, guide = F)
} else {
plt <- plt +
scale_fill_brewer(palette = "Set1", guide = F) +
scale_color_brewer(palette = "Set1", guide = F)
}
if (fix_yaxis) {
plt <- plt +
coord_cartesian(xlim = c(max(c(plot_start,
min(tib_gather$start))) / 1e6,
min(c(plot_end,
max(tib_gather$end))) / 1e6),
ylim = ylimits)
} else {
plt <- plt +
coord_cartesian(xlim = c(max(c(plot_start,
min(tib_gather$start))) / 1e6,
min(c(plot_end,
max(tib_gather$end))) / 1e6))
}
plot(plt)
}# Correlation plot only
metadata <- padamid_metadata_combined_act
tib_cor <- tib %>%
correlate(method = "pearson") %>%
gather(key, value, -term) %>%
mutate(value = replace_na(value, 1)) %>%
mutate(n1 = str_remove(term, "_Ki67"),
n2 = str_remove(key, "_Ki67"),
match1 = match(term, metadata$sample),
match2 = match(key, metadata$sample),
cell1 = metadata$cell[match1],
cell2 = metadata$cell[match2],
condition1 = metadata$condition[match1],
condition2 = metadata$condition[match2]) %>%
drop_na()
# Plot
plt <- tib_cor %>%
ggplot(aes(x = n1, y = n2, fill = value)) +
geom_tile() +
xlab("") + ylab("") +
scale_fill_gradientn(colors = colorRampPalette(rev(brewer.pal(n = 7,
name = "RdYlBu")))(100),
limits = c(min(-1, tib_cor$value), 1),
name = "Pearson correlation") +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
# Compare cell lines with each other before and after ActD
tib_cor %>%
filter(condition1 == condition2,
cell1 != cell2) %>%
ggplot(aes(x = cell1, y = value, col = cell2)) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_point() +
facet_grid(. ~ condition1) +
scale_color_manual(values = colors_set2[1:3]) +
xlab("") +
ylab("Pearson correlation") +
theme_bw() +
theme(aspect.ratio = 2,
axis.text.x = element_text(angle = 90, hjust = 1))
Modest correlations.
I also want to make these plots for individual chromosomes. Scatterplots are too much, but simple pearson correlations should do the trick.
# Calculate pearson correlations per chromosome
CorrelationsPerChromosome <- function(input_tib, metadata, ylim = c(-1.9, 1.9)) {
# Gather data and calculate pearson correlations
tib <- input_tib %>%
drop_na() %>%
dplyr::select(seqnames, metadata$sample)
# Calculate mean score per chromosome
tib_summary <- tib %>%
gather(key, value, -seqnames) %>%
mutate(match = match(key, metadata$sample),
cell = metadata$cell[match],
condition = metadata$condition[match],
target = metadata$target[match],
seqnames = factor(seqnames, chromosomes)) %>%
mutate(size = seqlengths(gr_padamid_combined)[match(seqnames,
seqlevels(gr_padamid_combined))],
size = size / 1e6,
rdna = seqnames %in% paste0("chr", c(13:15, 21:22))) %>%
drop_na()
# Boxplot of all scores per chromosome
plt <- tib_summary %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = value, fill = condition)) +
geom_boxplot(outlier.shape = NA) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(cell ~ interaction(target, condition)) +
xlab("") +
ylab("Score") +
coord_cartesian(ylim = ylim) +
scale_fill_manual(values = colors_set1, guide = F) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
# Single point (mean score) per chromosome
plt <- tib_summary %>%
group_by(cell, condition, target, size, seqnames) %>%
summarise(mean = mean(value)) %>%
ungroup() %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = mean, col = condition)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(target ~ cell) +
xlab("") +
ylab("Mean score") +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
# Single point (mean score) per chromosome
plt <- tib_summary %>%
group_by(cell, condition, target, size, seqnames, rdna) %>%
summarise(mean = mean(value)) %>%
ungroup() %>%
spread(condition, mean) %>%
mutate(`50ng` = `50ng` - DMSO) %>%
dplyr::select(-DMSO) %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames)),
cell = factor(cell, levels = levels(metadata$cell))) %>%
ggplot(aes(x = seqnames, y = `50ng`, col = cell, shape = rdna)) +
geom_point(size = 2) +
facet_grid(. ~ target) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
xlab("") +
ylab("Difference with DMSO") +
scale_color_manual(values = colors_set2) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
# Mean score vs chromosome size
plt <- tib_summary %>%
group_by(cell, condition, target, size, seqnames) %>%
summarise(mean = mean(value)) %>%
ggplot(aes(x = size, y = mean, col = condition)) +
geom_point() +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(cell ~ interaction(target, condition)) +
xlab("Chromosome size (Mb)") +
ylab("Mean score") +
coord_cartesian(xlim = c(0, 250)) +
scale_color_manual(values = colors_set1, guide = F) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
tib_cor <- tib %>%
group_by(seqnames) %>%
nest() %>%
mutate(cordf = map(data, correlate)) %>%
unnest(cordf) %>%
dplyr::select(-data) %>%
gather(key, value, -seqnames, -term) %>%
drop_na() %>%
mutate(value = unlist(value)) %>%
mutate(n1 = str_remove(term, "_Ki67"),
n2 = str_remove(key, "_Ki67"),
match1 = match(term, metadata$sample),
match2 = match(key, metadata$sample),
cell1 = metadata$cell[match1],
cell2 = metadata$cell[match2],
condition1 = metadata$condition[match1],
condition2 = metadata$condition[match2],
seqnames = factor(seqnames, chromosomes)) %>%
mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames,
seqlevels(gr_padamid_replicates))],
size = size / 1e6) %>%
filter(cell1 == cell2) %>%
drop_na() %>%
ungroup()
# Plot
plt <- tib_cor %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = value)) +
geom_point() +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(. ~ cell1) +
xlab("") +
ylab("Pearson correlation") +
coord_cartesian(ylim = c(min(0, tib_cor$value), 1)) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
plt <- tib_cor %>%
ggplot(aes(x = size, y = value)) +
geom_point() +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(. ~ cell1) +
xlab("Chromosome size (Mb)") +
ylab("Pearson correlation") +
coord_cartesian(xlim = c(0, 250),
ylim = c(min(0, tib_cor$value), 1)) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
plot(plt)
}
# Correlations per chromosome
#CorrelationsPerChromosome(tib_padamid_replicates, padamid_metadata_replicates_act)
CorrelationsPerChromosome(tib_padamid_combined, padamid_metadata_combined_act)
2. Enrichment at centromeres and telomeres
These analyses are based on simple counts. I want to do a similar analysis based on the normalized scores themselves.
# Get similar chromosomes for the two objects and filter samples
tib_padamid_act <- tib_padamid_combined %>%
dplyr::select(seqnames, start, end, padamid_metadata_combined_act$sample) %>%
filter(seqnames %in% chromosomes)
centromeres <- centromeres[seqnames(centromeres) %in% chromosomes]
# Also, prepare telomeres
telomeres <- c(GRanges(seqnames = chromosomes,
ranges = IRanges(start = 1, end = 1)),
GRanges(seqnames = chromosomes,
ranges = IRanges(start = seqlengths(gr_padamid_combined)[chromosomes],
end = seqlengths(gr_padamid_combined)[chromosomes])))
seqinfo(telomeres) <- seqinfo(gr_padamid_combined)
telomeres <- sort(telomeres)
# Add distance to centromeres and telomeres to the GRanges object - in Mb
dis <- distanceToNearest(as(tib_padamid_act, "GRanges"), centromeres)
tib_padamid_act$distance_to_centromere <- mcols(dis)$distance / 1e6
dis <- distanceToNearest(as(tib_padamid_act, "GRanges"), telomeres)
tib_padamid_act$distance_to_telomere <- mcols(dis)$distance / 1e6
# Process
tib <- tib_padamid_act %>%
dplyr::select(-start, -end) %>%
mutate(rdna = as.logical(seqnames %in% paste0("chr", c(13, 14, 15, 21, 22)))) %>%
gather(key, value, -contains("distance"), -rdna, -seqnames) %>%
separate(key, c("cell", "experiment", "condition", "target"), remove = F) %>%
mutate(match = match(key, padamid_metadata_combined_act$sample),
condition = padamid_metadata_combined_act$condition[match],
cell = padamid_metadata_combined_act$cell[match],
dis_group_centromere = as.numeric(cut(distance_to_centromere,
breaks = seq(-1, max(distance_to_centromere) + 5,
by = 1))) - 1,
dis_group_telomere = as.numeric(cut(distance_to_telomere,
breaks = seq(-1, max(distance_to_telomere) + 5,
by = 1))) - 2,
dis_group_telomere = ifelse(dis_group_telomere < 0, 0, dis_group_telomere)) %>%
dplyr::select(-contains("distance")) %>%
gather(class, class_value, contains("dis_group")) %>%
mutate(class = str_remove(class, "dis_group_"))
# Plot as boxplots
# First, calculate boxplots to have more control over ylimits
tib_summary <- tib %>%
group_by(cell, condition, class, class_value) %>%
drop_na() %>%
summarise(ymin = quantile(value, 0.05),
lower = quantile(value, 0.25),
middle = quantile(value, 0.5),
upper = quantile(value, 0.75),
ymax = quantile(value, 0.95)) %>%
ungroup() %>%
drop_na() %>%
arrange(cell, condition) %>%
mutate(sample = paste0(cell, "_", condition),
sample = factor(sample, levels = unique(sample)))
tib_summary %>%
filter(class_value <= 30) %>%
ggplot(aes(x = class_value, ymin = ymin, lower = lower, middle = middle,
upper = upper, ymax = ymax, group = class_value, y = middle)) +
geom_boxplot(stat="identity", outlier.shape = NA, fill = "lightgrey", col = "black") +
geom_line(aes(y = middle, group = NULL), col = "red") +
geom_hline(yintercept = 0, col = "blue", linetype = "dashed") +
facet_grid(class ~ sample, scales = "free_y") +
xlab("Distance to feature (Mb)") +
ylab("pA-DamID (log2)") +
#coord_cartesian(xlim = c(0, 30)) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
# Is the enrichment specific for rDNA chromosomes, or is chromosome size more
# important?
tib_summary <- tib %>%
filter(class_value < 2) %>%
mutate(condition = droplevels(condition)) %>%
group_by(cell, condition, seqnames, class) %>%
drop_na() %>%
summarise(mean = mean(value)) %>%
ungroup() %>%
mutate(seqnames = factor(seqnames, chromosomes),
rdna = seqnames %in% paste0("chr", c(13, 14, 15, 21, 22))) %>%
drop_na() %>%
mutate(size = seqlengths(gr_padamid_replicates)[match(seqnames,
seqlevels(gr_padamid_replicates))],
size = size / 1e6,
sample = paste0(cell, "_", condition),
sample = factor(sample, levels = unique(sample)))
tib_summary %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = mean, col = condition, shape = rdna)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(class ~ cell) +
xlab("") +
ylab("Ki67 score near centromeres") +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
# Add chromosome size and make a scatter plot
tib_summary %>%
ggplot(aes(x = size, y = mean, col = condition, shape = rdna)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
xlab("Chromosome size (Mb)") +
ylab("Ki67 score near centromeres (<2Mb)") +
scale_color_manual(values = colors_set3) +
facet_grid(class ~ cell) +
coord_cartesian(xlim = c(0, 250)) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, hjust = 1))
# Difference plot
tib_summary %>%
dplyr::select(-sample) %>%
spread(condition, mean) %>%
mutate(`50ng` = `50ng` - DMSO) %>%
dplyr::select(-DMSO) %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = size, y = `50ng`, shape = rdna, col = rdna)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(class ~ cell) +
xlab("") +
ylab("Ki67 score difference with DMSO near centromeres") +
coord_cartesian(xlim = c(0, 250)) +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
# Centromere only
tib_summary %>%
filter(class == "centromere") %>%
dplyr::select(-sample) %>%
spread(condition, mean) %>%
mutate(`50ng` = `50ng` - DMSO) %>%
dplyr::select(-DMSO) %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = size, y = `50ng`, shape = rdna, col = rdna)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(. ~ cell) +
xlab("") +
ylab("Ki67 score difference with DMSO near centromeres") +
coord_cartesian(xlim = c(0, 250)) +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
# Centromere only - ordered by chromosome size
tib_summary %>%
filter(class == "centromere") %>%
dplyr::select(-sample) %>%
spread(condition, mean) %>%
mutate(`50ng` = `50ng` - DMSO) %>%
dplyr::select(-DMSO) %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = `50ng`, shape = rdna, col = rdna)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(. ~ cell) +
xlab("") +
ylab("Ki67 score difference with DMSO near centromeres") +
#coord_cartesian(xlim = c(0, 250)) +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
# Centromere only - actual enrichments
tib_summary %>%
filter(class == "centromere") %>%
arrange(-size) %>%
mutate(seqnames = factor(seqnames, levels = unique(seqnames))) %>%
ggplot(aes(x = seqnames, y = mean, col = condition, shape = rdna)) +
geom_point(size = 2.5) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
facet_grid(. ~ cell) +
xlab("") +
ylab("Ki67 score near centromeres") +
scale_color_manual(values = colors_set1) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
This mostly shows the expected result: t=0h is enriched near telomeres while the other time points are more enriched at/near centromeres.
3. Example tracks
I want to make some clear example tracks of the differences in behaviour in the cell lines.
# Plot example tracks
# 1) RPE
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined,
exp_names = c("RPE_ActD_DMSO_Ki67",
"RPE_ActD_50ng_Ki67",
"RPE_ActD_DMSO_LMNB1",
"RPE_ActD_50ng_LMNB1",
"RPE_ActD_DMSO_H3K27me3",
"RPE_ActD_50ng_H3K27me3",
"RPE_ActD_DMSO_H3K9me3",
"RPE_ActD_50ng_H3K9me3"),
centromeres,
color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
RPE_ActD_50ng_Ki67 = "ActD",
RPE_ActD_DMSO_LMNB1 = "DMSO",
RPE_ActD_50ng_LMNB1 = "ActD",
RPE_ActD_DMSO_H3K27me3 = "DMSO",
RPE_ActD_50ng_H3K27me3 = "ActD",
RPE_ActD_DMSO_H3K9me3 = "DMSO",
RPE_ActD_50ng_H3K9me3 = "ActD"),
facet_groups = c(RPE_ActD_DMSO_Ki67 = "Ki67",
RPE_ActD_50ng_Ki67 = "Ki67",
RPE_ActD_DMSO_LMNB1 = "LMNB1",
RPE_ActD_50ng_LMNB1 = "LMNB1",
RPE_ActD_DMSO_H3K27me3 = "H3K27me3",
RPE_ActD_50ng_H3K27me3 = "H3K27me3",
RPE_ActD_DMSO_H3K9me3 = "H3K9me3",
RPE_ActD_50ng_H3K9me3 = "H3K9me3"),
smooth = 9, plot_chr = "chr2", fix_yaxis = F,
plot_start = 94e6, plot_end = 170e6,
color_list = colors_set1[c(1:5, 7)])## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
# 2) HCT116
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined,
exp_names = c("HCT116_ActD_DMSO_Ki67",
"HCT116_ActD_50ng_Ki67",
"HCT116_ActD_DMSO_LMNB1",
"HCT116_ActD_50ng_LMNB1",
"HCT116_ActD_DMSO_H3K27me3",
"HCT116_ActD_50ng_H3K27me3",
"HCT116_ActD_DMSO_H3K9me3",
"HCT116_ActD_50ng_H3K9me3"),
centromeres,
color_groups = c(HCT116_ActD_DMSO_Ki67 = "DMSO",
HCT116_ActD_50ng_Ki67 = "ActD",
HCT116_ActD_DMSO_LMNB1 = "DMSO",
HCT116_ActD_50ng_LMNB1 = "ActD",
HCT116_ActD_DMSO_H3K27me3 = "DMSO",
HCT116_ActD_50ng_H3K27me3 = "ActD",
HCT116_ActD_DMSO_H3K9me3 = "DMSO",
HCT116_ActD_50ng_H3K9me3 = "ActD"),
facet_groups = c(HCT116_ActD_DMSO_Ki67 = "Ki67",
HCT116_ActD_50ng_Ki67 = "Ki67",
HCT116_ActD_DMSO_LMNB1 = "LMNB1",
HCT116_ActD_50ng_LMNB1 = "LMNB1",
HCT116_ActD_DMSO_H3K27me3 = "H3K27me3",
HCT116_ActD_50ng_H3K27me3 = "H3K27me3",
HCT116_ActD_DMSO_H3K9me3 = "H3K9me3",
HCT116_ActD_50ng_H3K9me3 = "H3K9me3"),
smooth = 9, plot_chr = "chr2", fix_yaxis = F,
plot_start = 94e6, plot_end = 170e6,
color_list = colors_set1[c(1:5, 7)])## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
# 3) K562
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined,
exp_names = c("K562_ActD_DMSO_Ki67",
"K562_ActD_50ng_Ki67",
"K562_ActD_DMSO_LMNB1",
"K562_ActD_50ng_LMNB1",
"K562_ActD_DMSO_H3K27me3",
"K562_ActD_50ng_H3K27me3",
"K562_ActD_DMSO_H3K9me3",
"K562_ActD_50ng_H3K9me3"),
centromeres,
color_groups = c(K562_ActD_DMSO_Ki67 = "DMSO",
K562_ActD_50ng_Ki67 = "ActD",
K562_ActD_DMSO_LMNB1 = "DMSO",
K562_ActD_50ng_LMNB1 = "ActD",
K562_ActD_DMSO_H3K27me3 = "DMSO",
K562_ActD_50ng_H3K27me3 = "ActD",
K562_ActD_DMSO_H3K9me3 = "DMSO",
K562_ActD_50ng_H3K9me3 = "ActD"),
facet_groups = c(K562_ActD_DMSO_Ki67 = "Ki67",
K562_ActD_50ng_Ki67 = "Ki67",
K562_ActD_DMSO_LMNB1 = "LMNB1",
K562_ActD_50ng_LMNB1 = "LMNB1",
K562_ActD_DMSO_H3K27me3 = "H3K27me3",
K562_ActD_50ng_H3K27me3 = "H3K27me3",
K562_ActD_DMSO_H3K9me3 = "H3K9me3",
K562_ActD_50ng_H3K9me3 = "H3K9me3"),
smooth = 9, plot_chr = "chr2", fix_yaxis = F,
plot_start = 94e6, plot_end = 170e6,
color_list = colors_set1[c(1:5, 7)])## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
# 4) Effect of actinomycin D on Ki67 patterns
PlotDataTracksLight(input_tib = tib_padamid_combined,
exp_names = c("RPE_ActD_DMSO_Ki67",
"RPE_ActD_50ng_Ki67",
"HCT116_ActD_DMSO_Ki67",
"HCT116_ActD_50ng_Ki67",
"K562_ActD_DMSO_Ki67",
"K562_ActD_50ng_Ki67"),
centromeres,
color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
RPE_ActD_50ng_Ki67 = "ActD",
HCT116_ActD_DMSO_Ki67 = "DMSO",
HCT116_ActD_50ng_Ki67 = "ActD",
K562_ActD_DMSO_Ki67 = "DMSO",
K562_ActD_50ng_Ki67 = "ActD"),
smooth = 9, plot_chr = "chr2", fix_yaxis = T,
color_list = colors_set1[c(1:2)])## Warning: Removed 268 rows containing missing values (geom_rect).
PlotDataTracksLight(input_tib = tib_padamid_combined,
exp_names = c("RPE_ActD_DMSO_Ki67",
"RPE_ActD_50ng_Ki67",
"HCT116_ActD_DMSO_Ki67",
"HCT116_ActD_50ng_Ki67",
"K562_ActD_DMSO_Ki67",
"K562_ActD_50ng_Ki67"),
centromeres,
color_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
RPE_ActD_50ng_Ki67 = "ActD",
HCT116_ActD_DMSO_Ki67 = "DMSO",
HCT116_ActD_50ng_Ki67 = "ActD",
K562_ActD_DMSO_Ki67 = "DMSO",
K562_ActD_50ng_Ki67 = "ActD"),
smooth = 9, plot_chr = "chr22", fix_yaxis = T,
color_list = colors_set1[c(1:2)])## Warning: Removed 1472 rows containing missing values (geom_rect).
# 5) Combined tracks of cell types
# 3) K562
PlotDataTracksLightFaceted(input_tib = tib_padamid_combined,
exp_names = c("RPE_ActD_DMSO_Ki67",
"HCT116_ActD_DMSO_Ki67",
"K562_ActD_DMSO_Ki67",
"RPE_ActD_50ng_Ki67",
"HCT116_ActD_50ng_Ki67",
"K562_ActD_50ng_Ki67"),
centromeres,
color_groups = c(RPE_ActD_DMSO_Ki67 = "RPE",
HCT116_ActD_DMSO_Ki67 = "HCT116",
K562_ActD_DMSO_Ki67 = "K562",
RPE_ActD_50ng_Ki67 = "RPE",
HCT116_ActD_50ng_Ki67 = "HCT116",
K562_ActD_50ng_Ki67 = "K562"),
facet_groups = c(RPE_ActD_DMSO_Ki67 = "DMSO",
HCT116_ActD_DMSO_Ki67 = "DMSO",
K562_ActD_DMSO_Ki67 = "DMSO",
RPE_ActD_50ng_Ki67 = "ActD",
HCT116_ActD_50ng_Ki67 = "ActD",
K562_ActD_50ng_Ki67 = "ActD"),
smooth = 25, plot_chr = "chr2", fix_yaxis = F,
plot_start = 94e6, plot_end = 170e6,
color_list = colors_set2[c(1:3)])## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
4. Difference tracks
Let’s also make some difference tracks for IGV visualization.
MakeDifferenceTracks <- function(input_gr, metadata, basename, bigwig_dir,
target = "Ki67") {
# Function to create difference tracks of a certain experiment, using the
# basename to substract the values
# Get the sample names
samples <- metadata$sample
base_cell <- metadata$cell[metadata$sample == basename]
base_condition <- metadata$condition[metadata$sample == basename]
# Get the mcols for the granges
gr <- input_gr
mcols(gr) <- mcols(gr)[, samples]
# Prepare output directory
dir.create(bigwig_dir, showWarnings = F)
# Loop over the samples
for (s in samples) {
# Not for the basename
if (s == basename) next
# Not for different cell types
s_cell <- metadata$cell[metadata$sample == s]
s_condition <- metadata$condition[metadata$sample == s]
if (s_cell != base_cell) next
# Get the difference
tmp <- gr
mcols(tmp) <- data.frame(score = mcols(gr)[, s] - mcols(gr)[, basename])
tmp <- tmp[! is.na(tmp$score)]
# Save file
f_name <- paste0(s_cell, "_", s_condition, "_", target, "_vs_",
base_cell, "_", base_condition, "_", target, ".bw")
export.bw(tmp, file.path(bigwig_dir, f_name))
}
}
# Save tracks
bigwig_dir <- file.path(output_dir, "bigwig")
MakeDifferenceTracks(gr_padamid_combined, padamid_metadata_combined_act, "HCT116_ActD_DMSO_Ki67",
bigwig_dir)
MakeDifferenceTracks(gr_padamid_combined, padamid_metadata_combined_act, "K562_ActD_DMSO_Ki67",
bigwig_dir)5. Comparing the responses of Ki67 and K562 cells
The changes in Ki67 patterns upon actinomycin D treatment are very different between K562 and HCT116 cells. Why? Let’s correlate the differences with each other.
PlotScatter <- function(tib, n1, n2, color_by = NULL, identity = F) {
# Get tibble
tib <- tib %>%
dplyr::select(seqnames, matches(n1), matches(n2)) %>%
rename_all(~ c("seqnames", "n1", "n2"))
# Prepare color
if (! is.null(color_by)) {
tib <- tib %>%
add_column(color = color_by) %>%
drop_na()
alpha = 1
limits_color <- quantile(tib$color, c(0.001, 0.999), na.rm = T)
tib$color[tib$color < limits_color[1]] <- limits_color[1]
tib$color[tib$color > limits_color[2]] <- limits_color[2]
} else {
tib <- tib %>% drop_na()
tib$color = "1"
alpha = 0.02
}
# Plot
xlimits <- quantile(tib$n1, c(0.001, 0.999), na.rm = T) * 1.4
ylimits <- quantile(tib$n2, c(0.001, 0.999), na.rm = T) * 1.4
plt <- tib %>%
arrange(sample(1:nrow(.), size = nrow(.), replace = F)) %>%
ggplot(aes(x = n1, y = n2, color = color)) +
geom_point(size = 0.5, alpha = alpha) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
xlab(n1) +
ylab(n2) +
ggtitle(paste0("Spearman: ",
round(cor(tib$n1, tib$n2, use = "complete.obs",
method = "spearman"), 2))) +
coord_cartesian(xlim = xlimits, ylim = ylimits) +
theme_bw() +
theme(aspect.ratio = 1)
# Prepare color
if (! is.null(color_by)) {
plt <- plt +
scale_color_gradient2(low = "blue", mid = "grey", high = "red",
midpoint = 0)
} else {
plt <- plt +
scale_color_manual(values = "black", guide = F)
}
if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, col = "red", linetype = "dashed")
plot(plt)
}
# Get the differences
tib <- tib_padamid_combined %>%
mutate(RPE_difference = RPE_ActD_50ng_Ki67 - RPE_ActD_DMSO_Ki67,
HCT116_difference = HCT116_ActD_50ng_Ki67 - HCT116_ActD_DMSO_Ki67,
K562_difference = K562_ActD_50ng_Ki67 - K562_ActD_DMSO_Ki67,
cell_cycle_difference = RPE_6h_Ki67 - RPE_1h_Ki67)
# Plot the differences
set.seed(11)
PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_Ki67")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67", color_by = tib$HCT116_difference)
PlotScatter(tib, "HCT116_ActD_DMSO_Ki67", "K562_ActD_DMSO_Ki67", color_by = tib$K562_difference)
PlotScatter(tib, "RPE_difference", "HCT116_difference")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
PlotScatter(tib, "RPE_difference", "K562_difference")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
PlotScatter(tib, "HCT116_difference", "K562_difference")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
6. Comparison with RPE cell cycle time series
I have the hypothesis that actinomycin D “releases” Ki67 from the nucleoli and results in global Ki67 interactions. As Ki67 is much less limited to nucleoli early in interphase, the changes in Ki67 binding should thus correspond with the changes during the cell cycle. Possibly, this can be visualized by comparing the actinomycin D condition with the changes during the cell cycle.
# How do actinomycin D changes compare with cell cycle changes
PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_50ng_Ki67", identity = T)## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
PlotScatter(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_50ng_Ki67", identity = T,
color_by = tib$cell_cycle_difference)
PlotScatter(tib, "RPE_difference", "cell_cycle_difference")## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
## use `guide = "none"` instead.
Not really.
7. Effect of actinomycin D on H3K27me3, H3K9me3 and LMNB1
Show the effects on Ki67, H3K27me3, H3K9me3 and LMNB1. With individual scatter plots and a combined small plot.
# All scatterplots
tib_cor <- tibble()
cells <- c("RPE", "HCT116", "K562")
targets <- c("Ki67", "LMNB1", "H3K27me3", "H3K9me3")
for (cell in cells) {
for (target in targets) {
# Get the data
tib_tmp <- tib_padamid_combined %>%
dplyr::select(1:3,
matches("Act") &
matches(cell) &
matches(target)) %>%
rename_at(4:5, ~ c("DMSO", "ActD")) %>%
drop_na()
# Scatterplot
PlotScatterBinned(tib_tmp, "DMSO", "ActD", count_range = c(0, 900))
# Determine pearson correlation
cor_tmp <- cor(tib_tmp$DMSO, tib_tmp$ActD)
# Add this to tib_cor object
tib_cor <- bind_rows(tib_cor,
tibble(cell = cell,
target = target,
cor = cor_tmp))
}
}## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Add factors
tib_cor <- tib_cor %>%
mutate(cell = factor(cell, levels = cells),
target = factor(target, levels = targets))
# Combined plot
tib_cor %>%
ggplot(aes(x = target, y = cor, col = cell)) +
geom_point(size = 3) +
xlab("") +
ylab("Pearson correlation") +
scale_color_manual(values = colors_set2) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
theme_bw() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
And, where is the Ki67 binding, which chromatin marks?
# Ki67 patterns after nucleolar release
tib_cor <- tibble()
cells <- c("RPE", "HCT116", "K562")
targets <- c("LMNB1", "H3K27me3", "H3K9me3")
treatments <- c("DMSO", "50ng")
for (cell in cells) {
for (target in targets) {
for (treatment in treatments) {
# Get the data
tib_tmp <- tib_padamid_combined %>%
dplyr::select(1:3,
matches("Act") &
matches(cell) &
matches(treatment) &
matches(paste0("Ki67|", target))) %>%
rename_at(4:5, function(x) str_replace(x, ".*_", paste0(treatment, "_"))) %>%
drop_na()
# Scatterplot
PlotScatterBinned(tib_tmp,
paste0(treatment, "_", target),
paste0(treatment, "_", "Ki67"),
count_range = c(0, 750))
# Determine pearson correlation
cor_tmp <- cor(tib_tmp %>% pull(paste0(treatment, "_", target)),
tib_tmp %>% pull(paste0(treatment, "_", "Ki67")))
# Add this to tib_cor object
tib_cor <- bind_rows(tib_cor,
tibble(cell = cell,
target = target,
treatment = treatment,
cor = cor_tmp))
}
}
}## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Add factors
tib_cor <- tib_cor %>%
mutate(cell = factor(cell, levels = cells),
target = factor(target, levels = targets),
treatment = factor(treatment, levels = treatments))
# Combined plot
tib_cor %>%
ggplot(aes(x = target, y = cor, col = treatment, shape = treatment)) +
geom_point(size = 3) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
xlab("") +
ylab("Pearson correlation") +
scale_color_manual(values = colors_set1) +
facet_grid(. ~ cell) +
#scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
theme_bw() +
theme(aspect.ratio = 2,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
8. Detailed comparison of Ki67 changes versus LaminB1 and histones
# Calculate differences (smoothed)
tib <- tib_padamid_combined %>%
mutate(RPE_Ki67_diff = RPE_ActD_50ng_Ki67 - RPE_ActD_DMSO_Ki67,
RPE_LMNB1_diff = RPE_ActD_50ng_LMNB1 - RPE_ActD_DMSO_LMNB1,
RPE_H3K27me3_diff = RPE_ActD_50ng_H3K27me3 - RPE_ActD_DMSO_H3K27me3,
RPE_H3K9me3_diff = RPE_ActD_50ng_H3K9me3 - RPE_ActD_DMSO_H3K9me3,
HCT116_Ki67_diff = HCT116_ActD_50ng_Ki67 - HCT116_ActD_DMSO_Ki67,
HCT116_LMNB1_diff = HCT116_ActD_50ng_LMNB1 - HCT116_ActD_DMSO_LMNB1,
HCT116_H3K27me3_diff = HCT116_ActD_50ng_H3K27me3 - HCT116_ActD_DMSO_H3K27me3,
HCT116_H3K9me3_diff = HCT116_ActD_50ng_H3K9me3 - HCT116_ActD_DMSO_H3K9me3,
K562_Ki67_diff = K562_ActD_50ng_Ki67 - K562_ActD_DMSO_Ki67,
K562_LMNB1_diff = K562_ActD_50ng_LMNB1 - K562_ActD_DMSO_LMNB1,
K562_H3K27me3_diff = K562_ActD_50ng_H3K27me3 - K562_ActD_DMSO_H3K27me3,
K562_H3K9me3_diff = K562_ActD_50ng_H3K9me3 - K562_ActD_DMSO_H3K9me3) %>%
group_by(seqnames) %>%
mutate_at(vars(contains("diff")), function(x) runmean(x, k = 3)) %>%
ungroup()
# Binned scatter plots - LaminB1
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_DMSO_LMNB1",
n2 = "RPE_ActD_50ng_LMNB1",
color_by = tib$RPE_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "HCT116_ActD_DMSO_LMNB1",
n2 = "HCT116_ActD_50ng_LMNB1",
color_by = tib$HCT116_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "K562_ActD_DMSO_LMNB1",
n2 = "K562_ActD_50ng_LMNB1",
color_by = tib$K562_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Binned scatter plots - H3K27me3
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_DMSO_H3K27me3",
n2 = "RPE_ActD_50ng_H3K27me3",
color_by = tib$RPE_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "HCT116_ActD_DMSO_H3K27me3",
n2 = "HCT116_ActD_50ng_H3K27me3",
color_by = tib$HCT116_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "K562_ActD_DMSO_H3K27me3",
n2 = "K562_ActD_50ng_H3K27me3",
color_by = tib$K562_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Binned scatter plots - H3K9me3
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_DMSO_H3K9me3",
n2 = "RPE_ActD_50ng_H3K9me3",
color_by = tib$RPE_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "HCT116_ActD_DMSO_H3K9me3",
n2 = "HCT116_ActD_50ng_H3K9me3",
color_by = tib$HCT116_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib_padamid_combined, smooth = 3,
n1 = "K562_ActD_DMSO_H3K9me3",
n2 = "K562_ActD_50ng_H3K9me3",
color_by = tib$K562_Ki67_diff,
ylimits_col = c(-2.8, 2.8),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_Ki67_diff",
n2 = "HCT116_Ki67_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_Ki67_diff",
n2 = "K562_Ki67_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_Ki67_diff",
n2 = "K562_Ki67_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Ki67 vs LaminB1
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_Ki67_diff",
n2 = "RPE_LMNB1_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_Ki67_diff",
n2 = "HCT116_LMNB1_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_Ki67_diff",
n2 = "K562_LMNB1_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Ki67 vs H3K27me3
# RPE
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_Ki67_diff",
n2 = "RPE_H3K27me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_DMSO_Ki67",
n2 = "RPE_ActD_DMSO_LMNB1",
color_by = tib$RPE_ActD_DMSO_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_50ng_Ki67",
n2 = "RPE_ActD_50ng_LMNB1",
color_by = tib$RPE_ActD_50ng_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# HCT116
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_Ki67_diff",
n2 = "HCT116_H3K27me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_ActD_DMSO_Ki67",
n2 = "HCT116_ActD_DMSO_LMNB1",
color_by = tib$HCT116_ActD_DMSO_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_ActD_50ng_Ki67",
n2 = "HCT116_ActD_50ng_LMNB1",
color_by = tib$HCT116_ActD_50ng_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# K562
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_Ki67_diff",
n2 = "K562_H3K27me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_ActD_DMSO_Ki67",
n2 = "K562_ActD_DMSO_LMNB1",
color_by = tib$K562_ActD_DMSO_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_ActD_50ng_Ki67",
n2 = "K562_ActD_50ng_LMNB1",
color_by = tib$K562_ActD_50ng_H3K27me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# Ki67 vs H3K9me3
# RPE
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_Ki67_diff",
n2 = "RPE_H3K9me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_DMSO_Ki67",
n2 = "RPE_ActD_DMSO_LMNB1",
color_by = tib$RPE_ActD_DMSO_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "RPE_ActD_50ng_Ki67",
n2 = "RPE_ActD_50ng_LMNB1",
color_by = tib$RPE_ActD_50ng_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# HCT116
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_Ki67_diff",
n2 = "HCT116_H3K9me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_ActD_DMSO_Ki67",
n2 = "HCT116_ActD_DMSO_LMNB1",
color_by = tib$HCT116_ActD_DMSO_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "HCT116_ActD_50ng_Ki67",
n2 = "HCT116_ActD_50ng_LMNB1",
color_by = tib$HCT116_ActD_50ng_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# K562
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_Ki67_diff",
n2 = "K562_H3K9me3_diff", count_range = c(0, 1200),
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_ActD_DMSO_Ki67",
n2 = "K562_ActD_DMSO_LMNB1",
color_by = tib$K562_ActD_DMSO_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, smooth = 3,
n1 = "K562_ActD_50ng_Ki67",
n2 = "K562_ActD_50ng_LMNB1",
color_by = tib$K562_ActD_50ng_H3K9me3,
remove_outliers = T)## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
tib_cor <- tib %>%
drop_na() %>%
# Smoothing of 3 bins?
dplyr::summarise(RPE_LMNB1 = cor(RPE_Ki67_diff,
RPE_LMNB1_diff),
RPE_H3K27me3 = cor(RPE_Ki67_diff,
RPE_H3K27me3_diff),
RPE_H3K9me3 = cor(RPE_Ki67_diff,
RPE_H3K9me3_diff),
HCT116_LMNB1 = cor(HCT116_Ki67_diff,
HCT116_LMNB1_diff),
HCT116_H3K27me3 = cor(HCT116_Ki67_diff,
HCT116_H3K27me3_diff),
HCT116_H3K9me3 = cor(HCT116_Ki67_diff,
HCT116_H3K9me3_diff),
K562_LMNB1 = cor(K562_Ki67_diff,
K562_LMNB1_diff),
K562_H3K27me3 = cor(K562_Ki67_diff,
K562_H3K27me3_diff),
K562_H3K9me3 = cor(K562_Ki67_diff,
K562_H3K9me3_diff)) %>%
gather(key, value) %>%
separate(key, c("cell", "target"), remove = F) %>%
mutate(cell = factor(cell, levels(padamid_metadata_combined$cell)),
target = factor(target, levels(padamid_metadata_combined$target)))
tib_cor %>%
ggplot(aes(x = cell, y = value, col = target)) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
geom_point(size = 3) +
xlab("") +
ylab("Pearson correlation between differences") +
scale_color_manual(values = colors_set1[2:4]) +
theme_bw() +
theme(aspect.ratio = 1.5,
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
# DMSO
PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1",
color_by = tib$RPE_ActD_DMSO_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "RPE_ActD_DMSO_Ki67", "RPE_ActD_DMSO_LMNB1",
color_by = tib$RPE_ActD_DMSO_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1",
color_by = tib$HCT116_ActD_DMSO_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_DMSO_Ki67", "HCT116_ActD_DMSO_LMNB1",
color_by = tib$HCT116_ActD_DMSO_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1",
color_by = tib$K562_ActD_DMSO_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_DMSO_Ki67", "K562_ActD_DMSO_LMNB1",
color_by = tib$K562_ActD_DMSO_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
# 50ng ActD
PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1",
color_by = tib$RPE_ActD_50ng_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "RPE_ActD_50ng_Ki67", "RPE_ActD_50ng_LMNB1",
color_by = tib$RPE_ActD_50ng_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1",
color_by = tib$HCT116_ActD_50ng_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "HCT116_ActD_50ng_Ki67", "HCT116_ActD_50ng_LMNB1",
color_by = tib$HCT116_ActD_50ng_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1")## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1",
color_by = tib$K562_ActD_50ng_H3K27me3,
ylimits_col = c(-2.4, 2.2))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
PlotScatterBinned(tib, "K562_ActD_50ng_Ki67", "K562_ActD_50ng_LMNB1",
color_by = tib$K562_ActD_50ng_H3K9me3,
ylimits_col = c(-1.8, 1.8))## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.
Conclusions
Session info
sessionInfo()## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] knitr_1.33 caTools_1.18.2 corrr_0.4.3
## [4] GGally_2.1.2 RColorBrewer_1.1-2 rtracklayer_1.50.0
## [7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1
## [10] S4Vectors_0.28.1 BiocGenerics_0.36.1 forcats_0.5.1
## [13] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
## [16] readr_2.0.0 tidyr_1.1.3 tibble_3.1.6
## [19] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 bitops_1.0-7
## [3] matrixStats_0.60.0 fs_1.5.0
## [5] lubridate_1.7.10 httr_1.4.2
## [7] tools_4.0.5 backports_1.2.1
## [9] bslib_0.2.5.1 utf8_1.2.2
## [11] R6_2.5.1 mgcv_1.8-36
## [13] DBI_1.1.1 colorspace_2.0-2
## [15] withr_2.4.2 tidyselect_1.1.1
## [17] compiler_4.0.5 cli_3.1.0
## [19] rvest_1.0.1 Biobase_2.50.0
## [21] xml2_1.3.2 DelayedArray_0.16.3
## [23] labeling_0.4.2 sass_0.4.0
## [25] scales_1.1.1 digest_0.6.28
## [27] Rsamtools_2.6.0 rmarkdown_2.11
## [29] XVector_0.30.0 pkgconfig_2.0.3
## [31] htmltools_0.5.1.1 MatrixGenerics_1.2.1
## [33] highr_0.9 dbplyr_2.1.1
## [35] rlang_0.4.12 readxl_1.3.1
## [37] rstudioapi_0.13 farver_2.1.0
## [39] jquerylib_0.1.4 generics_0.1.0
## [41] jsonlite_1.7.2 BiocParallel_1.24.1
## [43] RCurl_1.98-1.3 magrittr_2.0.1
## [45] GenomeInfoDbData_1.2.4 Matrix_1.3-4
## [47] Rcpp_1.0.7 munsell_0.5.0
## [49] fansi_0.5.0 lifecycle_1.0.1
## [51] stringi_1.7.3 yaml_2.2.1
## [53] SummarizedExperiment_1.20.0 zlibbioc_1.36.0
## [55] plyr_1.8.6 grid_4.0.5
## [57] crayon_1.4.2 lattice_0.20-44
## [59] splines_4.0.5 Biostrings_2.58.0
## [61] haven_2.4.1 hms_1.1.0
## [63] pillar_1.6.4 codetools_0.2-18
## [65] reprex_2.0.0 XML_3.99-0.6
## [67] glue_1.5.0 evaluate_0.14
## [69] modelr_0.1.8 vctrs_0.3.8
## [71] tzdb_0.1.2 cellranger_1.1.0
## [73] gtable_0.3.0 reshape_0.8.8
## [75] assertthat_0.2.1 xfun_0.24
## [77] broom_0.7.9 GenomicAlignments_1.26.0
## [79] ellipsis_0.3.2